Décrire le problème à analyser

Nous chercherons à déterminer lequel des deux modèles (SARIMA ou Prophet) fournit les prévisions les plus précises et robustes pour cette série temporelle. Le but est de mieux comprendre les cycles de popularité du terme “Football” et de choisir le modèle le mieux adapté.

Appeler les packages nécessaires

library(prophet)
library(tseries)
library(forecast)
library(TSstudio)
library(dplyr)
library(plotly)
library(lmtest)

library(dygraphs)
library(lubridate)
library(Amelia)

Importer la série et la transformer à un objet de type ts()

serie2 <- read.csv("C:\\Users\\imadh\\Downloads\\multiTimeline.csv",sep=",")
data_p <- read.csv("C:\\Users\\imadh\\Downloads\\multiTimeline.csv",sep=",")


missmap(data_p) 
summary(data_p)

##      Date             Recherche     
##  Length:251         Min.   : 14.00  
##  Class :character   1st Qu.: 29.00  
##  Mode  :character   Median : 39.00  
##                     Mean   : 45.25  
##                     3rd Qu.: 57.00  
##                     Max.   :100.00

serie=serie2
serie=serie[,2]
class(serie)

## [1] "integer"

serie <- ts(serie, frequency = 12)
ts_info(serie)

##  The serie series is a ts object with 1 variable and 251 observations
##  Frequency: 12 
##  Start time: 1 1 
##  End time: 21 11

Tracer la série

ts.plot(serie,col="purple")

Diviser la série en partie apprentissage et partie test

data=ts_split(serie,sample.out = 12)
train=data$train
test=data$test
ts_info(train)

##  The train series is a ts object with 1 variable and 239 observations
##  Frequency: 12 
##  Start time: 1 1 
##  End time: 20 11

ts_info(test)

##  The test series is a ts object with 1 variable and 12 observations
##  Frequency: 12 
##  Start time: 20 12 
##  End time: 21 11

Transformation de la série(stationnariser la sarie)

# On applique une différenciation d'ordre 1 sur les données d'entraînement train. Cette opération est courante pour rendre une série stationnaire, en éliminant la tendance.
serie1=diff(train,k=1,differences = 1)

#on applique une différenciation supplémentaire sur serie1 avec un décalage de 12
serie2=diff(serie1,12,differences = 1)
ts.plot(serie1)
ts.plot(serie2)

Tester la stationnarité de la série transformée

adf.test(serie2)

## Warning in adf.test(serie2): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie2
## Dickey-Fuller = -7.6322, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

kpss.test(serie2)

## Warning in kpss.test(serie2): p-value greater than printed p-value

## 
##  KPSS Test for Level Stationarity
## 
## data:  serie2
## KPSS Level = 0.029421, Truncation lag parameter = 4, p-value = 0.1

Sarima

Le modèle SARIMA, qui signifie Seasonal Autoregressive Integrated Moving Average (moyenne mobile autoregressive intégrée saisonnière), est un modèle de prévision de séries temporelles polyvalent et largement utilisé. C’est une extension du modèle ARIMA non saisonnier, conçu pour traiter des données avec des motifs saisonniers. SARIMA capte à la fois les dépendances à court terme et à long terme dans les données, ce qui en fait un outil robuste pour la prévision. Il combine les concepts de modèles autorégressifs (AR), intégrés (I), et de moyenne mobile (MA) avec des composantes saisonnières.

• Les Composantes de SARIMA

Pour comprendre SARIMA, examinons ses différentes composantes :

– Composante Saisonnière

Le “S” dans SARIMA représente la saison. La saisonnalité fait référence à des motifs répétitifs dans les données, pouvant être quotidiens, mensuels, annuels, ou à toute autre période régulière. Identifier et modéliser la composante saisonnière est l’un des points forts de SARIMA.

– Composante Autorégressive (AR)

La partie AR dans SARIMA représente la composante autorégressive, qui modélise la relation entre la valeur actuelle et ses valeurs passées. Elle capture l’autocorrélation des données, c’est-à-dire comment les données sont corrélées avec elles-mêmes au fil du temps.

– Composante Intégrée (I)

La partie I dans SARIMA indique le différenciation, qui transforme les données non stationnaires en données stationnaires. La stationnarité est cruciale pour la modélisation des séries temporelles. La composante intégrée mesure combien de différences sont nécessaires pour rendre la série stationnaire.

– Composante Moyenne Mobile (MA)

La partie MA dans SARIMA représente la composante de moyenne mobile, qui modélise la dépendance entre la donnée actuelle et les erreurs de prédiction passées. Elle aide à capter le bruit à court terme dans les données.

• Différenciation Saisonnière

Avant de plonger dans SARIMA, il est essentiel de comprendre la différenciation saisonnière. La différenciation saisonnière est le processus de soustraction des données de série temporelle par un décalage égal à la saisonnalité. Cela permet de supprimer la composante saisonnière et de rendre les données stationnaires, ce qui facilite la modélisation. La différenciation saisonnière est souvent notée par “D” dans SARIMA.

• Notation SARIMA

La notation SARIMA est la suivante :

où : - AR(p) : Composante autorégressive d’ordre p. - MA(q) : Composante de moyenne mobile d’ordre q. - I(d) : Composante intégrée d’ordre d. - AR saisonnier(P) : Composante autorégressive saisonnière d’ordre P. - MA(Q) : Composante de moyenne mobile saisonnière d’ordre Q. - I(D) : Composante intégrée saisonnière d’ordre D. - s : Période saisonnière.

• Représentation mathématique de SARIMA

La représentation mathématique du modèle SARIMA est la suivante :

où : - $y_t$ est la série temporelle observée au temps $t$, - $B$ est l’opérateur de retard, représentant le décalage temporel (c’est-à-dire $B y_t = y_{t-1}$), - $\varphi_1$ est le coefficient autorégressif non saisonnier, - $\Phi_1$ est le coefficient autorégressif saisonnier, - $\theta_1$ est le coefficient de moyenne mobile non saisonnier, - $\Theta_1$ est le coefficient de moyenne mobile saisonnier, - $s$ est la période saisonnière, - $\epsilon_t$ est le terme d’erreur aléatoire au temps $t$.

• Composantes de l’équation

– Composante Autorégressive (AR) : La composante autorégressive non saisonnière, représentée par $(1 - \varphi_1 B)$, capture la relation entre l’observation actuelle et certaines valeurs retardées (précédentes dans la série temporelle). Le terme $B$ représente l’opérateur de décalage, qui déplace la série temporelle en arrière par une certaine période.

– Composante Autorégressive Saisonnière (SAR) : La composante autorégressive saisonnière est représentée par $(1 - \Phi_1 B^s)$. Cette composante capture la relation entre l’observation actuelle et certaines valeurs retardées à des intervalles saisonniers.

– Composante de Différenciation Non Saisonnière : La composante de différenciation non saisonnière est représentée par $(1 - B)$, où $d$ est l’ordre de la différenciation non saisonnière. Cette composante est utilisée pour rendre la série stationnaire en la différenciant un certain nombre de fois.

– Composante de Différenciation Saisonnière : La composante de différenciation saisonnière est représentée par $(1 - B^s)$, où $D$ est l’ordre de la différenciation saisonnière. Cette composante est utilisée pour rendre la série stationnaire en la différenciant à des intervalles saisonniers.

– Composante de Moyenne Mobile (MA) : La composante de moyenne mobile est représentée par $(1 + \theta_1 B)$. Cette composante capture la relation entre l’observation actuelle et les erreurs résiduelles issues d’un modèle de moyenne mobile appliqué aux observations retardées.

– Composante de Moyenne Mobile Saisonnière (SMA) : La composante de moyenne mobile saisonnière est représentée par $(1 + \Theta_1 B^s)$. Elle capture la relation entre l’observation actuelle et les erreurs résiduelles issues d’un modèle de moyenne mobile appliqué aux observations retardées à des intervalles saisonniers.

– Terme d’Erreur : Le terme d’erreur, $\epsilon_t$, représente le bruit aléatoire ou la variation inexpliquée dans la série temporelle.

# Set up the ACF and PACF plots to examine series autocorrelations

par(mfrow=c(1,2))
acf(serie2,lag.max = 64 ,main = "ACF of Training Series")
pacf(serie2,lag.max = 64,main = "PACF of Training Series")

# Define the range of parameters for model selection
p=q=P=Q=0:2
grid=expand.grid(p,q,P,Q)
names(grid)=c("p","q","P","Q")
grid$d = 1  # Difference for non-seasonal component
grid$D = 1  # Difference for seasonal component
head(grid)

##   p q P Q d D
## 1 0 0 0 0 1 1
## 2 1 0 0 0 1 1
## 3 2 0 0 0 1 1
## 4 0 1 0 0 1 1
## 5 1 1 0 0 1 1
## 6 2 1 0 0 1 1

# Add a column to sum parameter orders for each model
grid$k=rowSums(grid)
# Filter the grid for models with parameter sum k <= 6
grid=grid%>%filter(k<=6)
head(grid)

##   p q P Q d D k
## 1 0 0 0 0 1 1 2
## 2 1 0 0 0 1 1 3
## 3 2 0 0 0 1 1 4
## 4 0 1 0 0 1 1 3
## 5 1 1 0 0 1 1 4
## 6 2 1 0 0 1 1 5

dim(grid)

## [1] 50  7

# Initialize AIC column
grid$aic=0

# Loop through each combination of parameters in the grid
for(i in 1:nrow(grid)){
 md=arima(train,order=c(grid$p[i],1,grid$q[i]),
 seasonal=list(order=c(grid$P[i],1,grid$Q[i])))
 grid$aic[i]=md$aic
}
# Arrange models by AIC in ascending order and select the model with the lowest AIC
grid=arrange(grid,aic)
# Fit the model with the best parameters
i=1
md=arima(train,order=c(grid$p[i],1,grid$q[i]),
 seasonal=list(order=c(grid$P[i],1,grid$Q[i])))
## afficher le modèle estimé
md

## 
## Call:
## arima(x = train, order = c(grid$p[i], 1, grid$q[i]), seasonal = list(order = c(grid$P[i], 
##     1, grid$Q[i])))
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.6011  -1.0000  -0.6967
## s.e.  0.0537   0.0175   0.0499
## 
## sigma^2 estimated as 34.69:  log likelihood = -728.5,  aic = 1465

coeftest(md)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.601145   0.053682  11.198 < 2.2e-16 ***
## ma1  -0.999999   0.017487 -57.184 < 2.2e-16 ***
## sma1 -0.696667   0.049917 -13.957 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

checkresiduals(md) #

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(0,1,1)[12]
## Q* = 13.118, df = 21, p-value = 0.9044
## 
## Model df: 3.   Total lags used: 24

qqnorm(md$residuals) #diagramme quantil-quantil des résidus
jarque.bera.test(md$residuals) # test de normalité de jarque bera

## 
##  Jarque Bera Test
## 
## data:  md$residuals
## X-squared = 539.91, df = 2, p-value < 2.2e-16

library(lmtest)
coeftest(md)# tester la sig des coefs

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1   0.601145   0.053682  11.198 < 2.2e-16 ***
## ma1  -0.999999   0.017487 -57.184 < 2.2e-16 ***
## sma1 -0.696667   0.049917 -13.957 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fc=forecast(md,h=12) #prédire les 12 observations futures
accuracy(fc,test) # calculer les indicateurs statistiques

##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.2400237 5.727543 3.656115 -2.916995 8.823205 0.6780541
## Test set      2.6322461 4.964846 3.299041  6.404850 7.541021 0.6118320
##                    ACF1 Theil's U
## Training set 0.02625275        NA
## Test set     0.47830448 0.5013221

# comparer l'estimation du modèle avec la vraie série
test_forecast(serie,forecast.obj = fc,test=test) 

Estimer le modèle avec la série entière

#Ajuste un modèle ARIMA sur la série entière avec des paramètres fixes
i=1
sarima_forecast=arima(serie,order=c(grid$p[i],1,grid$q[i]),
 seasonal=list(order=c(grid$P[i],1,grid$Q[i])))
sarima_forecast

## 
## Call:
## arima(x = serie, order = c(grid$p[i], 1, grid$q[i]), seasonal = list(order = c(grid$P[i], 
##     1, grid$Q[i])))
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.6031  -1.0000  -0.7114
## s.e.  0.0522   0.0165   0.0468
## 
## sigma^2 estimated as 33.52:  log likelihood = -763.02,  aic = 1534.04

on a obtenit un modèle avec un meilleur ajustement

Prédire les 12 observations futures

fc=forecast(sarima_forecast,12)
plot(fc)
plot_forecast(fc)

Estimer le modèle à l’aide de auto.arima()

# max.order représente: p+q+P+Q
md2 <- auto.arima(train,
 max.order = 5,
 D = 1,
 d = 1)
md2

## Series: train 
## ARIMA(0,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ma1      ma2     sma1
##       -0.3307  -0.4028  -0.6747
## s.e.   0.0656   0.0754   0.0489
## 
## sigma^2 = 37.56:  log likelihood = -732.98
## AIC=1473.96   AICc=1474.15   BIC=1487.65

fc1=forecast(md2,h=12)
test_forecast(serie,forecast.obj = fc1,test=test)

plot_forecast(fc1)

checkresiduals(md2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 18.57, df = 21, p-value = 0.6127
## 
## Model df: 3.   Total lags used: 24

qqnorm(md2$residuals)
jarque.bera.test(md2$residuals)

## 
##  Jarque Bera Test
## 
## data:  md2$residuals
## X-squared = 244.31, df = 2, p-value < 2.2e-16

Prophet Facebook

Le modele Facebook Prophet repose sur un modèle décomposable, c’est-à-dire qu’il sépare le comportement d’une série temporelle en plusieurs composantes : la tendance, la saisonnalité et les jours fériés. Cela permet de créer des modèles prédictifs plus précis en prenant en compte ces différentes contraintes. Il est souvent considéré comme supérieur au modèle ARIMA, car il aide à ajuster et affiner les paramètres d’entrée.

Équation mathématique du modèle Prophet

L’équation du modèle Prophet est la suivante :

où :

- $y(t)$ fait référence à la prévision.

- $g(t)$ représente la tendance.

- $s(t)$ représente la saisonnalité.

- $h(t)$ fait référence aux jours fériés pour la prévision.

- $e(t)$ représente l’erreur de prévision.

Termes importants utilisés dans le modèle Facebook Prophet

• Tendance

  La tendance désigne un changement dans le développement de la série temporelle, que ce soit dans une direction croissante ou décroissante. Mathématiquement, elle est représentée par :

où :

- $C$ est la capacité de saturation.

- $k$ est le taux de croissance.

- $m$ est le paramètre de décalage.

• Seasonality

  La Seasonality fait référence à une caractéristique d’une série temporelle qui se répète à un moment spécifique ou lors d’une saison particulière, ce qui modifie la tendance. La saisonnalité peut être annuelle, hebdomadaire, ou même quotidienne, selon la nature des données.

• Holidays

  Les Holidays sont des périodes qui peuvent avoir un impact significatif sur l’activité d’une entreprise. Ils peuvent entraîner des variations dans les profits ou les pertes, selon le secteur d’activité. L’inclusion des jours fériés dans le modèle permet d’ajuster la prévision pour ces événements spéciaux qui perturbent les tendances normales.

# Création du dataframe pour Prophet
prophet_data <- data_p %>%
  rename(ds = Date, y = Recherche)

# Vérifier les premières lignes pour s'assurer que ds est bien au format Date
head(prophet_data)

##           ds  y
## 1 2004-01-01 34
## 2 2004-02-01 30
## 3 2004-03-01 29
## 4 2004-04-01 30
## 5 2004-05-01 31
## 6 2004-06-01 37

# Ajuster le modèle Prophet
prophet_model <- prophet::prophet(prophet_data, weekly.seasonality = FALSE, daily.seasonality = FALSE)

# Créer le dataframe pour les prédictions futures
future <- prophet::make_future_dataframe(prophet_model, periods = 12, freq = "month")

# Forecast
forecast_prophet <- predict(prophet_model, future)

# Afficher les prévisions
tail(forecast_prophet[c('ds', 'yhat', 'yhat_lower', 'yhat_upper')])

##             ds     yhat yhat_lower yhat_upper
## 258 2025-06-01 34.77448   24.85727   43.94813
## 259 2025-07-01 34.61974   26.12873   43.88786
## 260 2025-08-01 51.79199   42.93789   61.00854
## 261 2025-09-01 83.43805   73.89502   92.19512
## 262 2025-10-01 72.26260   62.71667   81.60312
## 263 2025-11-01 71.01135   60.80925   80.52490

# Plot les prévisions
dyplot.prophet(prophet_model, forecast_prophet)

Reset Zoom

Actual

Predicted

10

20

30

40

50

60

70

80

90

100

2010

2020

# Plot the trend
plot(forecast_prophet$ds, forecast_prophet$trend, type = "l", main = "Trend Component", col = "blue")

# Plot the seasonal effects
plot(forecast_prophet$ds, forecast_prophet$yearly, type = "l", main = "Yearly Seasonality", col = "red")

1. MSE (Mean Squared Error):

• SARIMA : 13.87

• Prophet : 52.93

• Le MSE est une mesure de l’erreur quadratique moyenne entre les valeurs observées et les prévisions. Un MSE plus faible est préférable. Ici, SARIMA a un MSE beaucoup plus faible, ce qui suggère que le modèle SARIMA a des erreurs de prédiction plus petites en moyenne comparé à Prophet.

2. RMSE (Root Mean Squared Error):

• SARIMA : 3.72

• Prophet : 7.28

• Le RMSE est la racine carrée du MSE et est également utilisé pour mesurer la magnitude des erreurs. Une valeur plus faible est meilleure. SARIMA a un RMSE inférieur, ce qui indique que les prévisions de SARIMA sont plus précises en termes d’écart par rapport aux valeurs réelles que celles de Prophet.

3. MAE (Mean Absolute Error):

• SARIMA : 2.54

• Prophet : 5.08

• Le MAE mesure l’erreur absolue moyenne, c’est-à-dire la différence moyenne entre les valeurs observées et les prévisions. Moins l’erreur est élevée, mieux c’est. Ici encore, SARIMA surpasse Prophet, car il produit une erreur absolue plus faible en moyenne.

4. MPE (Mean Percentage Error):

• SARIMA : -5.40%

• Prophet : 2.67%

• Le MPE mesure l’erreur en pourcentage, où une valeur proche de 0% est idéale. Une valeur négative indique une tendance à sous-estimer, tandis qu’une valeur positive indique une tendance à surestimer les prévisions. SARIMA sous-estime légèrement la série temporelle, tandis que Prophet a tendance à surestimer.

5. MAPE (Mean Absolute Percentage Error):

• SARIMA : 5.89%

• Prophet : 12.27%

• Le MAPE est un indicateur de l’erreur en pourcentage absolu. Un MAPE plus bas indique une meilleure performance du modèle. SARIMA a un MAPE beaucoup plus bas que Prophet, ce qui indique que SARIMA génère des prévisions plus précises (en pourcentage de l’erreur absolue) que Prophet.

6. MASE (Mean Absolute Scaled Error):

• SARIMA : 2.54

• Prophet : 0.50

• Le MASE est une métrique qui permet de comparer un modèle par rapport à un modèle de référence (par exemple, un modèle naïf). Un MASE inférieur à 1 indique que le modèle est meilleur que le modèle naïf. Prophet a un MASE inférieur à 1, ce qui signifie qu’il est plus performant que le modèle de base, tandis que SARIMA, avec un MASE plus élevé, n’est pas aussi bon comparé à un modèle de référence naïf.

Conclusion:

• SARIMA semble mieux fonctionner pour cette série temporelle particulière, car il a un MSE, un RMSE, un MAE, un MAPE plus bas que Prophet, ce qui signifie que ses prévisions sont plus précises et qu’il fait moins d’erreurs.

• Prophet, bien que plus simple à utiliser et plus adapté à la prise en compte des effets saisonniers et des jours fériés, semble avoir de plus grandes erreurs en termes absolus et relatifs pour cette série temporelle, ce qui peut indiquer qu’il n’est pas aussi bien adapté à ce jeu de données particulier.